Dalton Lecture: How Far Can We Go In Distributed Hydrological Modelling?

This paper considers distributed hydrological models in hydrology as an expression of a pragmatic realism. Some of the problems of distributed modelling are discussed including the problem of nonlinearity, the problem of scale, the problem of equifinality, the problem of uniqueness and the problem of uncertainty. A structure for the application of distributed modelling is suggested based on an uncertain or fuzzy landscape space to model space mapping. This is suggested as the basis for an Alternative Blueprint for distributed modelling in the form of an application methodology. This Alternative Blueprint is scientific in that it allows for the formulation of testable hypotheses. It focusses attention on the prior evaluation of models in terms of physical realism and on the value of data in model rejection. Finally, some unresolved questions are outlined that distributed modelling must address in the future together with a vision for distributed modelling as a means of learning about places. Realism in the face of adversity It is almost 30 years since I wrote my first distributed hydrological model for my PhD thesis, following the Freeze and Harlan (1969) blueprint but using finite element methods. My thesis (Beven, 1975) contained an application of the model to the small East Twin catchment in the UK, the same catchment that had been studied in the field by Weyman (1970). The model represented a catchment as a number of variable width, slope following, hillslope segments, each represented by a 2D (vertical and downslope directions) solution of the Richards equation (Figure 1). Computer limitations meant that only a coarse finite element mesh could be used and that even then, on the computers available, it proved difficult to perform simulations that took less computer time than real time simulated. The modelling results were never published. They were simply not good enough. The model did not reproduce the stream discharges, it did not reproduce the measured water table levels, it did not reproduce the observed heterogeneity of inputs into the stream from the hillslopes (Figure 2). It was far easier at the time to publish the results of hypothetical simulations (Beven, 1977). The ideas in what follows are essentially a distillation of those early experiences and of thinking hard about how to do distributed modelling in some sense “properly” since then. The limitations of that PhD study were in part because of the crudeness of the representation given the computer resources available at the time (the model itself actually existed as two boxes of computer cards). Just as in numerical weather forecasting, the accuracy of numerical algorithms for solving the partial differential equations and the feasible discretisation of the flow domains has improved dramatically since 1975. However, just as in numerical weather forecasting, there remain limits to the detail that can be represented and there remains a problem of representing or parameterising sub-grid scale processes. As computer power improves further into the future, the feasible discretisation will become finer but the problem of sub-grid parameterisation does not go away. The form of that parameterisation might become simpler at finer scale but there is then the problem of knowing what the actual values of parameters for all the different spatial elements might be (Beven, 1989, 1996b, 2000a). There is then an interesting question as to how far such models, with their necessary approximations of processes and parameters at the element scale, can represent reality. An analysis of this question reveals a number of issues. These will be summarised here as the problems of nonlinearity; of scale; of uniqueness; of equifinality; and of uncertainty. The aim is, as ever, a “realistic” representation of the hydrology of a catchment that will be useful in making predictions in situations that have not yet occurred or where measurements have not yet been made. Indeed, one argument for the use of distributed modelling in hydrology has always been that they might be more “realistic” than simpler models that are calibrated to historical data in a curve-fitting exercise, with no guarantee, therefore, that they might do well in simulating responses in other periods or other conditions (e.g. Beven and O’Connell, 1982, Beven, 1985). That argument continues to be used in discussions of the problems of parameter estimation (e.g. Smith et al., 1994; De Marsily, 1994; Beven et al., 2001). What then does “realism” mean in the context of distributed hydrological modelling. At the risk of making a gross generalisation I would suggest that most practising environmental scientists have, as a working philosophy, a pragmatic or heuristic realism; that the quantities that we deal with exist independently of our perceptions and empirical studies of them, that this extends even to quantities that are not (yet) observable, and that further work will move the science towards a more realistic description of the world. Again at the risk of generalising I would also suggest that most practising environmental scientists do not worry too much about the theory laden nature of their studies, (subsuming any such worries within the general framework of the critical rationalist stance that things will get better as studies progress). As has been pointed out many times this theory laden-ness applies very much to experimental work, but it applies even more pointedly to modelling work where theory must condition model results very strongly. This pragmatic realism is a “natural” philosophy in part because, as environmental scientists we are often dealing with phenomena that are close to our day-to-day perceptions of the world. At a fundamental level I do a lot of computer modelling but I think of it as representing real water. If I try to predict pollutant transport I think of it as trying to represent a real pollutant. Environmental chemists measure the characteristics of real solutions and so on. What I am calling pragmatic realism naturally combines elements of objectivism, actualism, empiricism, idealism, instrumentalism, Bayesianism, relativism and hermeneutics; of multiple working hypotheses, falsification, and critical rationalism (but allowing adjustment of auxiliary conditions); of confirmation and limits of validity; of methodologies of research programmes while maintaining an open mind to paradigm shifts; and of the use of “scientific method” within the context of the politics of grant awarding programmes and the sociology of the laboratory. Refined and represented in terms of ideals rather than practice, it probably comes closest to the transcendental realism of Bhaskar (1989, Collier, 1994). However, in hydrology, at least, the practice often appears to have more in common with the entertaining relativism of Feyerabend (1991), not least because theories are applied to systems that are open which, as Cartwright (1999) has recently pointed out even makes the application of the equation force=mass*acceleration difficult to verify or apply in practice in many situations. Hydrologists also know only too well the difficulties of verifying or applying the mass and energy balance equations in open systems (Beven, 2001b, d). This does not, of course, mean that such principles or laws should not be applied in practice, only that we should be careful about the limitations of their domain of validity (as indeed are engineers in the application of the force equation). It is in the critical rationalist idea that the description of reality will continue to improve that many of the problems of environmental modelling have been buried for a long time. This apparent progress is clearly the case in many areas of environmental science such as weather forecasting and numerical models of the ocean. It is not nearly so clear in distributed hydrological modelling even though many people feel that, by analogy, it should be. This analogy is critically misguided, for some of the reasons that will be explored in the sections that follow. It has led to a continuing but totally unjustified determinism in many applications of distributed modelling and a lack of recognition of how far we can go in distributed hydrological modelling in the face of these adverse problems. The problem of nonlinearity The problem of nonlinearity is at the heart of many of the problems faced in the application of distributed modelling concepts in hydrology, despite the fact that for many years “linear” models, such as the unit hydrograph and more recent linear transfer functions, have been shown to work well (see, for example, Beven 2001a), particularly in larger catchments (but see Goodrich et al., 1995, for a counter-example in a semi-arid environment where channel transmission losses result in greater apparent nonlinearity with increasing catchment size). In fact, this apparent linearity is a de facto artefact of the analysis. It applies only to the relationship between some “effective” rainfall inputs and river discharge (and sometimes only to the “storm runoff” component of discharge). It does not apply to the relationship between rainfall inputs and river discharge that is known to be a nonlinear function of antecedent conditions, rainfall volume, and the (interacting) surface and subsurface processes of runoff generation. Hydrological systems are nonlinear and the implications of this nonlinearity should be taken into account in the formulation and application of distributed models. This we do attempt to do, of course. All distributed models have nonlinear functional relationships included in their local element scale process descriptions of surface and subsurface runoff generation, whether they are based on Richards equation or the SCS curve number. We have not been so good at taking account of some of the other implications of dealing with nonlinear dynamical systems, however. These include, critically, the fact that nonlinear equations do not average simply and that the extremes of any distribution of responses in a nonlinear system may be important in controlling the observed responses. Crudely interpreted in hydrological terms, this means local subgrid-scale nonlinear descriptions, such as Richards equation, should not be used at the model element scale (let alone at the GCM grid scale) where the heterogeneity of local parameter variations is expected to be important (Beven, 1989). The local heterogeneities will mean that the element scale averaged equations must be different from the local scale descriptions; that using mean local scale parameter values will not give the correct results, especially where there are coupled surface and subsurface flows (Binley et al., 1989); and that the extremes of the local responses (infiltration rates, preferential flows, areas of first saturation) will be important. This suggests, for example, that the use of pedotransfer functions to estimate a set of average soil parameters at the element scale of a distributed hydrological model should not be expected to give accurate results. Note: this follows purely from considerations of nonlinear mathematics, even if Richards equations is acceptable as a description of the local flow processes (which could also be debated, e.g. Beven and Germann, 1982). These implications are well known, so why have they been ignored for so long in distributed modelling in hydrology? Is it simply because there is no “physically based” theory to put in the place of Richards equation, since alternative sub-grid paramterisations seem too “conceptual” in nature? The recent work by Reggiani et al. (1998, 1999, 2000) is an attempt to formulate equations directly at the subcatchment or flow element scale directly in terms of mass, energy and momentum equations but has not solved the problem of parameterising the integrated exchanges between elements in heterogeneous flow domains. There are other implications of nonlinearity that are known to be important. Nonlinear systems are sensitive to their initial and boundary conditions. Unconstrained they will often exhibit chaotic behaviour. Initial and boundary conditions are poorly known in hydrology (see notably Stephenson and Freeze, 1974), as often are the observed values with which we compare model predictions, but fortunately the responses are necessarily constrained by mass and energy balances. It is these constraints that have allowed hydrological modellers to avoid worrying too much about the potential for chaos. Essentially, by maintaining approximately correct mass and energy balances our models cannot go too far wrong, especially after a bit of calibration of parameter values. That does not mean, however, that it is easy to get very good predictions (even allowing for observation error), especially for extreme events. This is reinforced by recent work in nonlinear dynamics looking at stochastically forced systems of simple equations. This work suggests that where there is even a slight error in the behaviour or attractor of a model of the system, the model will not be able to correctly reproduce the extremes of the distribution of the output variables either for short time scales or for integrated outputs over long (e.g. annual) time scales. If this is true for simple systems, does it imply that the same should be true for flood prediction and water yield predictions using (always slightly wrong) distributed models in hydrology? How can predictive capability be protected against these effects of nonlinearity? The problem of scale The problem of scale is inherently linked to that of nonlinearity. Scale issues in linear systems are only related to the problem of adequately assessing the inputs at different scales with available measurement techniques. As is well known by all hydrological modellers this is a problem even in the simple assessment of rainfalls over different sizes of catchment area, even before trying to make some assessment of the nature and heterogeneity of the surface and subsurface processes with the measurement techniques available. It is clear, for example, that we have kept the Richards equation approach as a subgrid scale parameterisation for so long because it is consistent with the measurement scales of soil physical measurements. Because we have no measurement techniques that directly give information at the element grid scales (say 10m to 1km in the case of distributed hydrological models to 5 to 100km in the case of NWP and GCM models) we have not developed the equivalent, scale consistent, process descriptions that would then implicitly take account of the effects of subgrid scale heterogeneity and nonlinearity. A recent comment by Blöschl (2001) has discussed the scale problem in hydrology. His analysis has much the same starting point as that of Beven (1995). He also recognises the need to identify the “dominant process controls” at different scales but comes to a totally different conclusion. Whereas Beven (1995) suggests that scaling theories will ultimately prove to be impossible and that is therefore necessary to recognise the scale dependence of model structures, Blöschl (2001) suggested that it is in resolving the scale problem that the real advances will be made in hydrological theorising and practice in the future. How do these two viewpoints bear on the application of distributed hydrological models? Let us assume for the moment that it might be possible to develop a scaling theory that would allow the definition of grid or element scale equations and parameter values on the basis of knowledge of the parameter values at smaller scales. Certainly some first attempts have been made to do so in subsurface flows (Dagan, 1986, and others) and surface runoff (Tayfur and Kavvas, 1998). Attempts are also being made to describe element scale processes in terms of more fundamental characteristics of the flow domain, such as depositional scenarios for sedimentary aquifers. This reveals the difference between hydrology and some other subject areas in this respect. In hydrology, the development of a scaling theory is not just a matter of the dynamics and organisation of the flow of the fluid itself. In surface and subsurface hillslope hydrology, the flow is always responding to the local pore scale or surface boundary conditions. The characteristics of the flow domain determine the flow velocities. Those characteristics must be represented as parameter values at some scale. Those parameter values must be estimated in some way. But the characteristics are impossible to determine everywhere, even for surface runoff if it occurs. For subsurface flow processes the characteristics are essentially unknowable with current measurement techniques. Thus they must be inferred in some way from either indirect or large scale measurements. In both cases, a theory of inference would be required. This would be the scaling theory but it is clear from this argument that any such theory would need to be supported by strong assumptions about the nature of the characteristics of the flow domain even if we felt secure about the nonlinearities of the flow process descriptions. The assumptions would not, however, be verifiable: in fact it is more likely that they would be made for mathematical tractability rather than physical realism and applied without being validated for a particular flow domain because, again, of the limitations of current measurement techniques. Thus, the problem of scale in distributed hydrological modelling does not arise because we do not know the principles involved. We do, if we think about it, understand a lot about the issues raised by nonlinearities of the processes, heterogeneities of the flow domains, limitations of measurement techniques, and the problem of knowing parameter values or structures everywhere. The principles are general and we have at least a qualitative understanding of their implications, but the difficulty comes in the fact that we are required to apply hydrological models in particular catchments, all with their own unique characteristics. The problem of uniqueness In the last 30 years of distributed hydrological modelling there has been an implicit underlying theme of developing a general theory of hydrological processes. It has been driven by the pragmatic realist philosophy outlined earlier. The idea that if we can get the description of the dynamics of the processes correct then parameter identification problems will become more tractable is still strongly held. However, in a recent paper, I have put forward an alternative view: that we should take much more account of the particular characteristics of particular catchment areas, i.e. the uniqueness of place (Beven, 2000a). It is useful in this respect to consider the case where we could define the “perfect” model description. In its equations, such a model would properly reflect all the effects of local heterogeneity on the flow dynamics and the nonlinearities associated with the coupling of different flow processes. Test simulations with such a model would show how it takes account of the redistribution of the inputs by a vegetation cover; the initiation of local overland flows, reinfiltration on heterogeneous surfaces, initiation and propagation of preferential flows etc. Such a model clearly has the potential to produce predictions that are accurate to within the limitations of measurement errors. However, such a model must still have some way of taking account of all the local heterogeneities of the flow domain in any application to a particular catchment. In short, even the perfect model has parameters that have to be estimated. Presumably, the perfect model will embody within it some expressions to relate the parameter values it requires to some measureable characteristics of the flow domain (indeed, the perfect model seems to require that a scaling theory is, in fact, feasible). This could be done in either a disaggregation or aggregation framework. A disaggregation framework would require making inferences from catchment scale measurements to smaller scale process parameters. This would be similar to the type of calibration exercise against catchment discharges that is often carried out today. It clearly leaves scope for multiple parameter sets being able to reproduce the catchment scale behaviour in a way that is consistent with the model dynamics. An aggregation process implies that information will be required on the heterogeneity of parameter values within the catchment area. We will not, however, be able to determine those parameters everywhere in a particular catchment area with its own unique characteristics, especially because the perfect model would tell us that it is the extremes of the distribution of characteristics that may be important in controlling storm runoff generation. It is always more difficult to estimate the extremes of a distribution than the first two moments. Thus a very large number of measurements would be required without any real guarantee that they might be spatially coherent. Since our current measurement techniques have severe limitations in assessing spatial variability then it would seem that the aggregation approach would also result in a large number of model parameter sets being consistent with the model dynamics in reproducing the large scale behaviour. Thus, even if we knew the structure of the perfect model, uniqueness of place leads to a very important identifiability problem. In the case of the perfect model, this could be considered as simply a problem of non-identifiability i.e. a unique (“optimal”) set of parameters would exist, if only we had the measurements available to be able to identify it. In practice, with limited measurements available there would most probably be a non-uniqueness problem i.e. that there appear to be several or many different optimal parameter sets but the measurements do not allow us to distinguish between them. However, we cannot normally assume that we are using such a perfect model structure. Thus, Beven (1993, 1996a,b) has suggested that it is better to approach the problem of uniqueness of place using a concept of equifinality of model structures and parameter sets. This choice of word is intended to indicate an explicit recognition that, given the limited measurements available in any application of a distributed hydrological model, it will not be possible to identify an “optimal” model. Rather, we should accept that there may be many different model structures and parameter sets that will be acceptable in simulating the available data. It is worth stressing in this that, even if we believed that we knew the perfect model structure, it would not be immune to the problem of equifinality in applications to particular catchments with their own unique characteristics. Limited measurements, and particularly the unknowability of the subsurface, will result in equifinality, even for the perfect model. There has been a commonly expressed hope that, in the future, remote sensing information would lead to the possibility of more robust estimates of spatially distributed parameter values for distributed hydrological modelling in applications to unique catchment areas. Pixel sizes for remote sensing are at the same scale, or even sometimes finer, than distributed model element scales and in many images we can easily detect visually spatial patterns that appear to be hydrologically significant (we can include here ground probing radar and cross-borehole tomography techniques that give some insight into the local nature of the subsurface flow domain). However, the potential for remote sensing to provide the information required would appear to be limited. The digital numbers stored by the sensor do not give direct estimates of the hydrogical variables or parameters required at the pixel scale. They require an interpretative model. Such a model will, itself, require parameter values to reflect the nature of the surface, the structure and state of the vegetation, the state of the atmosphere, etc. In fact, the digital numbers received by the user may already have been processed by an interpretative model to correct for atmospheric effects etc. in a way that may not reflect all the processes involved even if the interpretative model is physically “realistic”. The user may wish to leave such corrections to the imaging “experts”, but will then need to apply a further interpretative model for the hydrological purposes he/she has in mind. The resulting uncertainties may, at least sometimes, be very significant (see for example Franks et al., 1997), especially where the parameters of the interpretative model might also be expect to change over time, e.g. with vegetation growth or senescence. Thus, remote sensing information will also be subject to equifinality in interpretation and uncertainty in prediction. This will be compounded by the need to couple interpretative models for satellite or aircraft images which, except under unusual circumstances, give only information on near surface emissions, to models of the subsurface. However, it is worth repeating that it is often possible to see hydrologically significant patterns in some images. Thus, it should be expected that there is useful information on the distributed responses of particular hillslopes and catchments to be gained from remote sensing, but it will certainly not solve the problem of uniqueness. The problem of equifinality The recognition of equifinality arose out of Monte Carlo experiments in applying models with different parameter sets in simulating catchment scale discharges (Beven and Binley, 1992; Duan et al., 1992; Beven, 1993). It resulted in some interestingly different responses. The University of Arizona group response was that a better method for identifying the optimal parameter set was required leading to their development of the stochastic complex evolution methodology, as embodied in the UA-SCE algorithm. Other experiments in global optimisation have explored simulated annealing, genetic algorithms and Monte Carlo Markov Chain methodologies (e.g. Kuczera, 1997, Kuczera and Parent, 1999). A further recognition that the results of even a global optimisation depended strongly on the evaluation measure used has lead to the exploration of multi-objective optimisation techniques such as the Pareto optimal set methodology of Yapo et al. (1998) and Gupta et al. (1999), again from the Arizona group. The underlying aim, however, has still been to identify parameter sets that are in some sense optimal. The response of the Lancaster University group was different. They were prepared to reject the idea that an optimal model would ever be identifiable and develop the concept of equifinality in a more direct way. This lead to the Generalised Likelihood Uncertainty Estimation (GLUE) Methodology (Beven and Binley, 1992; Beven et al., 2000, Beven, 2001a). GLUE is an extension of the Generalised Sensitivity Analysis of Hornberger, Spear and Young (Hornberger and Spear, 1981; Spear et al., 1994) in which many different model parameter sets are chosen randomly, simulations run, and evaluation measures used to reject some models (model structure/parameter set combinations) as non-behavioural while all those considered as behavioural are retained in prediction. In GLUE the predictions of the behavioural models are weighted by a likelihood measure based on past performance to form a cumulative weighted distribution of any predicted variable of interest. Traditional statistical likelihood measures can be used in this framework, in which case the output prediction distributions can be considered as probabilities of prediction of the variable of interest. However, the methodology is general in that more general likelihood measures, including fuzzy measures, can be used in which case only conditional prediction limits or possibilities are estimated. Different likelihood measures can be combined using Bayes equation or a number of other methods (Beven et al., 2000; Beven, 2001a). There is one other implication of equifinality that is of particular importance in distributed modelling. Distributed models have the potential to use different parameter values for every different element in the spatial discretisation. In general this means that many hundreds or thousands of parameter values must be specified. Clearly it is not possible to optimise all these parameter values, they must be estimated on the basis of some other information, such as soil texture, vegetation type, surface cover etc. Values are available for different types of soil, vegetation etc in the literature. However, such values will themselves have been back-calculated or optimised against observations gathered in specific (unique) locations under particular sets of forcing conditions. One of the lessons from GLUE studies is that it is the parameter set that is important in giving a good fit to the observations. It is very rarely the case that the simulations are so sensitive to a particular parameter that only certain values of that parameter will give good simulations. More often a particular parameter value will give either good or bad simulations depending on the other parameter values in the set. Thus, bringing together different parameter values from different sources is no guarantee that, even if they were optimal in the situations where they were determined, they will give good results as a set in a new set of